L10: Geometric operations with rasters

Mosaicking, Resampling, Reprojecting

Bogdan G. Popescu

bogdan.popescu@johncabot.edu

John Cabot University

Mosicking Rasters

In this lectures, we will look at making changes to geometric component of rasters:

mosaicking
resampling
reprojecting

We will use packages like stars and units

Mosaicking: Elevation

In the next we few slides, we will use a Digitial Elevation Models (DEM) of the world, by mosaicking, subsetting, and reprojecting

We start putting all tiles together

Downloading DEM data

DEM stands for a digital elevation model

There are many data sources

The most helpful website for downloading DEM is: http://www.webgis.com/terr_world.html

Working with Elevation

Elevation

#Step1: Load the DEM files for every continent
#EUROPE
filepath<-"./data/elevation/w020n90/W020N90.DEM"
x1<-read_stars(filepath)
filepath<-"./data/elevation/w020n40/w020n40.DEM"
x2<-read_stars(filepath)
filepath<-"./data/elevation/e020n40/e020n40.DEM"
x3<-read_stars(filepath)
filepath<-"./data/elevation/e020n90/e020n90.DEM"
x4<-read_stars(filepath)
filepath<-"./data/elevation/w060n40/w060n40.DEM"
x5<-read_stars(filepath)
filepath<-"./data/elevation/w060n90/w060n90.DEM"
x6<-read_stars(filepath)

#ASIA
filepath<-"./data/elevation/e140n40/e140n40.DEM"
x7<-read_stars(filepath)
filepath<-"./data/elevation/e140n90/e140n90.DEM"
x8<-read_stars(filepath)
filepath<-"./data/elevation/e060n40/e060n40.DEM"
x9<-read_stars(filepath)
filepath<-"./data/elevation/e060n90/e060n90.DEM"
x10<-read_stars(filepath)
filepath<-"./data/elevation/e100n40/e100n40.DEM"
x11<-read_stars(filepath)
filepath<-"./data/elevation/e100n90/e100n90.DEM"
x12<-read_stars(filepath)

#AUSTRALIA
filepath<-"./data/elevation/e060s10/e060s10.DEM"
x13<-read_stars(filepath)
filepath<-"./data/elevation/e100s10/e100s10.DEM"
x14<-read_stars(filepath)
filepath<-"./data/elevation/e140s10/e140s10.DEM"
x15<-read_stars(filepath)
filepath<-"./data/elevation/e060s60/e060s60.DEM"
x16<-read_stars(filepath)
filepath<-"./data/elevation/e120s60/e120s60.DEM"
x17<-read_stars(filepath)

#NORTH AMERICA
filepath<-"./data/elevation/w180n90/w180n90.DEM"
x18<-read_stars(filepath)
filepath<-"./data/elevation/w140n90/w140n90.DEM"
x19<-read_stars(filepath)
filepath<-"./data/elevation/w100n90/w100n90.DEM"
x20<-read_stars(filepath)
filepath<-"./data/elevation/w180n40/w180n40.DEM"
x21<-read_stars(filepath)
filepath<-"./data/elevation/w140n40/w140n40.DEM"
x22<-read_stars(filepath)
filepath<-"./data/elevation/w100n40/w100n40.DEM"
x23<-read_stars(filepath)

#SOUTH AMERICA
filepath<-"./data/elevation/w180s10/w180s10.DEM"
x24<-read_stars(filepath)
filepath<-"./data/elevation/w140s10/w140s10.DEM"
x25<-read_stars(filepath)
filepath<-"./data/elevation/w100s10/w100s10.DEM"
x26<-read_stars(filepath)
filepath<-"./data/elevation/w180s60/w180s60.DEM"
x27<-read_stars(filepath)
filepath<-"./data/elevation/w120s60/w120s60.DEM"
x28<-read_stars(filepath)

#SOUTH AFRICA
filepath<-"./data/elevation/w060s10/w060s10.DEM"
x29<-read_stars(filepath)
filepath<-"./data/elevation/w020s10/w020s10.DEM"
x30<-read_stars(filepath)
filepath<-"./data/elevation/e020s10/e020s10.DEM"
x31<-read_stars(filepath)
filepath<-"./data/elevation/w060s60/w060s60.DEM"
x32<-read_stars(filepath)
filepath<-"./data/elevation/w000s60/w000s60.DEM"
x33<-read_stars(filepath)

Elevation

#Step2: Creating a mosaic
star_mosaic_eu <- st_mosaic(x1, x2, x3, x4, x5, x6)
star_mosaic_rus <- st_mosaic(x7, x8, x9, x10, x11, x12)
star_mosaic_aus <- st_mosaic(x13, x14, x15, x16, x17)
star_mosaic_namer <- st_mosaic(x18, x19, x20, x21, x22, x23)
star_mosaic_samer <- st_mosaic(x24, x25, x26, x27, x28)
star_mosaic_safr <- st_mosaic(x29, x30, x31, x32, x33)

#Step3: Creating a second mosaic
all<-st_mosaic(star_mosaic_eu, star_mosaic_rus, star_mosaic_aus, star_mosaic_namer, star_mosaic_samer, star_mosaic_safr)

#Step4: Higher dx/day values means lower resolution
elev<-st_warp(all, st_as_stars(st_bbox(all), dx = 0.4, dy = 0.4))

#Step5: Renaming the raster file
names(elev)<-"elev"

#Step6: Ensuring the the raster has the same crs
st_crs(elev)<-st_crs(borders)

fig<-ggplot()+
  geom_stars(data=elev)+
  scale_fill_gradientn(colours=c("chartreuse3","red"),   na.value = NA)+
  geom_sf(data=borders, linewidth = 0.01, fill = NA, color = "white", alpha=0.5)+
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

Digital Elevation Model (DEM)

Otherwise we only have on tile

Digital Elevation Model (DEM) in Europe

This is how we identify the coordinates of our frame.

library(purrr)

#Identifying coordinates of Northern points
norway<-subset(borders, SOVEREIGNT == "Norway")
nc_geom <- st_geometry(norway)
list_points_norway<-nc_geom[[1]]
list_points_norway<-flatten(list_points_norway)
list_points_norway2<-as.data.frame(do.call(rbind, list_points_norway))
names(list_points_norway2)<-c("lon_x", "lat_y")
max_lat_y_eu<-max(list_points_norway2$lat_y)

#Identifying coordinates of Southern points
greece<-subset(borders, SOVEREIGNT == "Greece")
nc_geom <- st_geometry(greece)
list_points_greece<-nc_geom[[1]]
list_points_greece<-flatten(list_points_greece)
list_points_greece2<-as.data.frame(do.call(rbind, list_points_greece))
names(list_points_greece2)<-c("lon_x", "lat_y")
min_lat_y_eu<-min(list_points_greece2$lat_y)

#Identifying coordinates of Eastern points
ukraine<-subset(borders, SOVEREIGNT == "Ukraine")
nc_geom <- st_geometry(ukraine)
list_points_ukraine<-nc_geom[[1]]
list_points_ukraine<-flatten(list_points_ukraine)
list_points_ukraine2<-as.data.frame(do.call(rbind, list_points_ukraine))
names(list_points_ukraine2)<-c("lon_x", "lat_y")
max_lon_x_eu<-max(list_points_ukraine2$lon_x)

#Identifying coordinates of Western points
portugal<-subset(borders, SOVEREIGNT == "Portugal")
nc_geom <- st_geometry(portugal)
list_points_portugal<-nc_geom[[1]]
list_points_portugal<-flatten(list_points_portugal)
list_points_portugal2<-as.data.frame(do.call(rbind, list_points_portugal))
names(list_points_portugal2)<-c("lon_x", "lat_y")
min_lon_x_eu<-min(list_points_portugal2$lon_x)

Digital Elevation Model (DEM) in Europe

#Step1: Higher dx/day values mean lower resolution
all2<-st_warp(all, st_as_stars(st_bbox(all), dx = 0.1, dy = 0.1))
#Step2: Renaming raster values
names(all2)<-"elev"
#Step3: Mapping
fig<-ggplot()+
  geom_stars(data=all2)+
  scale_fill_gradientn(colours=c("chartreuse3","red"),   na.value = NA)+
  geom_sf(data=borders, linewidth = 0.05, fill = NA, color = "white", alpha=0.5)+
  coord_sf(xlim = c(min_lon_x_eu-3, max_lon_x_eu+3), ylim = c(min_lat_y_eu-1, max_lat_y_eu-6), expand = FALSE)+
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

Digital Elevation Model (DEM) in Europe

Digital Elevation Model (DEM) in the US

We can do the same thing for the US

#Step1: Higher dx/day values mean lower resolution
all2<-st_warp(all, st_as_stars(st_bbox(all), dx = 0.1, dy = 0.1))
#Step2: Renaming raster values
names(all2)<-"elev"
#Step3: Identifying the coordinates for the map
min_lon_x_us<-(-130)
max_lon_x_us<-(-57)
min_lat_y_us<-(26)
max_lat_y_us<-(53)
#Step4: Mapping
fig<-ggplot()+
  geom_stars(data=all2)+
  scale_fill_gradientn(colours=c("chartreuse3","red"),   na.value = NA)+
  geom_sf(data=borders, linewidth = 0.05, fill = NA, color = "white", alpha=0.5)+
  coord_sf(xlim = c(min_lon_x_us-3, max_lon_x_us+3), ylim = c(min_lat_y_us-3, max_lat_y_us+3), expand = FALSE)+
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

Digital Elevation Model (DEM) in the US

We can do the same thing for the US:

Elevation 3D

#Step1: Turning the raster into points
elev2<-st_as_sf(elev, as_points = FALSE, merge = FALSE)
#Step2: Mapping
fig1<-ggplot() +
  geom_sf(data = elev2, aes(fill = elev), color=NA) +
  scale_fill_gradient(low = "chartreuse3", high = "red")+
  geom_sf(data=borders, linewidth = 0.01, fill = NA, color = "white", alpha=0.5)+
  theme_bw()
#Step3: Creating the movie
library("rayshader")
library("av")
plot_gg(fig1,multicore=TRUE,width=7,height=3,scale=350)    # Plot_gg de rayshader
render_snapshot(filename = "Elevation")
render_movie("./figures/movie_elev.mp4",frames = 720, fps=30,zoom=0.8,fov = 30)

Elevation 3D

Raster Resampling

Raster resampling is the process of transferring raster values from the original grid to a different grid.

Resampling may be necessary for:

alligning several input rasters that come from different sources to the same grid
Reducing the resolution of very detailed rasters, so that they are more convenient to work with

Raster Resampling

In case you have not noticed, we have already resampled the DEM raster.

#Step2: Creating a mosaic
star_mosaic_eu <- st_mosaic(x1, x2, x3, x4, x5, x6)
star_mosaic_rus <- st_mosaic(x7, x8, x9, x10, x11, x12)
star_mosaic_aus <- st_mosaic(x13, x14, x15, x16, x17)
star_mosaic_namer <- st_mosaic(x18, x19, x20, x21, x22, x23)
star_mosaic_samer <- st_mosaic(x24, x25, x26, x27, x28)
star_mosaic_safr <- st_mosaic(x29, x30, x31, x32, x33)

#Step3: Creating a second mosaic
all<-st_mosaic(star_mosaic_eu, star_mosaic_rus, star_mosaic_aus, star_mosaic_namer, star_mosaic_samer, star_mosaic_safr)

#Step4: Higher dx/day values mean lower resolution
elev<-st_warp(all, st_as_stars(st_bbox(all), dx = 0.4, dy=0.4))

#Step5: Renaming the raster file
names(elev)<-"elev"

#Step6: Ensuring the the raster has the same crs
st_crs(elev)<-st_crs(borders)

fig<-ggplot()+
  geom_stars(data=elev)+
  scale_fill_gradientn(colours=c("chartreuse3","red"),   na.value = NA)+
  geom_sf(data=borders, linewidth = 0.01, fill = NA, color = "white", alpha=0.5)+
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

Nearest Neighbor Resampling

There are a variety of other resampling methods

In the nearest neighbor sampling, the new raster pixels get the value from the nearest pixel of the original raster.

Some values may be lost.

# Nearest Neighbor Resampling {visibility=“uncounted”}

#Step4: Lower dx values. Lower dx values means lower resolution
grid = st_as_stars(st_bbox(all), dx = 0.4, dy = 0.4)
dem1 = st_warp(elev, grid, method = "near")

fig<-ggplot()+
  geom_stars(data=dem1)+
  scale_fill_gradientn(colours=c("chartreuse3","red"),   na.value = NA)+
  geom_sf(data=borders, linewidth = 0.01, fill = NA, color = "white", alpha=0.5)+
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
fig

Bilinear resampling

In bilinear resampling, each new raster cell gets a weighted average of four nearest cells from the input, rather than just one.

Bilinear resampling

Two technical things to note about st_warp when using a method other than method="near":

use_gdal=TRUE needs to be specified, otherwise the method argument is ignored
no_data_value needs to be set to a value outside of the valid range (such as no_data_value=-9999 for a DEM raster);

dem2 = st_warp(elev, grid, method = "bilinear", use_gdal = TRUE, no_data_value = -9999)

Average resampling

Another useful method is the average resampling method, where each new cell gets the weighted average of all overlapping input cells:

Average resampling

dem3 = st_warp(elev, grid, method = "average", use_gdal = TRUE, no_data_value = -9999)

Raster Reprojection

Raster reprojection is more complex than vector layer reprojection.

It requires both transforming pixel coordinates and a resampling step in order to “arrange” the transformed values back into a regular grid.

While other sources recommend st_warp for reprojection, this a lot of pressure on your computer’s memory.

Raster Reprojection

#Step1: Reprojecting the shapefile to EPSG:32636
borders_rep<-st_transform(borders, crs = 32636)
#Step2: Subsetting countries
italy<-subset(borders_rep, SOVEREIGNT=="Italy")
#Step3: Reducing the original resolution: Higher dx/day values mean lower resolution
elev<-st_warp(all, st_as_stars(st_bbox(all), dx = 0.4, dy=0.4))
#Step4: Turning the raster to a raster object 
dem_rep <- terra::rast(elev)
#Step5: Re-projecting the rasters
dem_rep<-terra::project(dem_rep, terra::crs(borders_rep))
#Step6: Turning the raster back to a stars object
dem_rep<-st_as_stars(dem_rep)

#The following section is to indetify the min max of Italy in ordet to map it.
#Step7: Making a box object around the polygon
bbox <- st_bbox(italy)

#Step8: Extracting min and max latitudes and longitudes
min_lat_y <- bbox["ymin"]
max_lat_y <- bbox["ymax"]
min_lon_x <- bbox["xmin"]
max_lon_x <- bbox["xmax"]

#Step9: Adding an error
error<-10000000
#Step10: Mapping
fig<-ggplot()+
  geom_stars(data=dem_rep)+
  scale_fill_gradientn(colours=c("chartreuse3","red"),   na.value = NA)+
  #geom_sf(data=italy, linewidth = 0.01, fill = NA, color = "black", alpha=0.5)
  geom_sf(data=italy, linewidth = 0.1, fill = NA)+
  geom_sf(data=borders_rep, linewidth = 0.1, fill = NA)+
  coord_sf(xlim = c(min_lon_x - error, max_lon_x + error), 
           ylim = c(min_lat_y - error, max_lat_y + error))

Raster Reprojection

fig

Raster Reprojection

Note that the coordinate units in the reprojected raster are no longer degrees, but meters.

We can check that by analyzing the output from st_crs(dem_rep)

st_crs(dem_rep)

Coordinate Reference System:
  User input: WGS 84 / UTM zone 36N 
  wkt:
PROJCRS["WGS 84 / UTM zone 36N",
    BASEGEOGCRS["WGS 84",
        ENSEMBLE["World Geodetic System 1984 ensemble",
            MEMBER["World Geodetic System 1984 (Transit)"],
            MEMBER["World Geodetic System 1984 (G730)"],
            MEMBER["World Geodetic System 1984 (G873)"],
            MEMBER["World Geodetic System 1984 (G1150)"],
            MEMBER["World Geodetic System 1984 (G1674)"],
            MEMBER["World Geodetic System 1984 (G1762)"],
            MEMBER["World Geodetic System 1984 (G2139)"],
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ENSEMBLEACCURACY[2.0]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["UTM zone 36N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",33,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Engineering survey, topographic mapping."],
        AREA["Between 30°E and 36°E, northern hemisphere between equator and 84°N, onshore and offshore. Belarus. Cyprus. Egypt. Ethiopia. Finland. Israel. Jordan. Kenya. Lebanon. Moldova. Norway. Russian Federation. Saudi Arabia. Sudan. Syria. Turkey. Uganda. Ukraine."],
        BBOX[0,30,84,36]],
    ID["EPSG",32636]]

Combining Raster and Vector Layers

Now, we will focus on elements related to:

Crop and mask a raster according to a vector layer
Switching from vector to raster representation, and vice versa
Calculate a raster of distances to nearest point
Extract raster values from locations defined by a vector laye

Masking and Cropping Rasters

Masking a raster measing turning pixels values outside of an area of interest—defined using a polygonal layer—into NA.

Masking and Cropping Rasters

#Step1: Only choosing Italy
italy<-subset(borders, SOVEREIGNT=="Italy")
#Step2: Reducing the resolution: Higher dx/day values mean lower resolution
elev<-st_warp(all, st_as_stars(st_bbox(all), dx = 0.08, dy = 0.08))
#Step3: Cropping
dem_italy<-st_crop(elev, italy)
#Step4: Renaming the raster file
names(dem_italy)<-"elev"

fig<-ggplot()+
  geom_stars(data=dem_italy)+
  scale_fill_gradientn(colours=c("chartreuse3","red"), na.value = NA)+
  geom_sf(data=italy, linewidth = 0.1, fill = NA)+
  theme_bw()
fig

Turning a Raster into a Vector

It is possible to turn rasters into vectors

We have already seen this in fact.

library(sf)
#Step0: Switched off spherical geometry (s2)
sf::sf_use_s2(FALSE)
#Step1: Only choosing Italy
italy<-subset(borders, SOVEREIGNT=="Italy")
#Step1: Ensuring the same CRS for vector and raster
st_crs(italy)<-st_crs(borders)
st_crs(elev)<-st_crs(borders)
#Step2: Cropping raster
dem_italy<-st_crop(elev, italy)
#Step3: Changing the name
names(dem_italy)<-"elev"
#Step3: Creating a grid from the raster
dem_italy2<-st_as_sf(dem_italy, as_points = FALSE)
#Step4: Changing the name of the raster
names(dem_italy2)[1]<-"elev"

Turning a Raster into a Vector

This fundamentally creates a grid

#Step7: Making a box object around the polygon
bbox <- st_bbox(italy)

#Step8: Extracting min and max latitudes and longitudes
min_lat_y <- bbox["ymin"]
max_lat_y <- bbox["ymax"]
min_lon_x <- bbox["xmin"]
max_lon_x <- bbox["xmax"]

error<-1
fig1<-ggplot()+
  geom_sf(data=dem_italy2, aes(fill=elev), linewidth=0)+
  geom_sf(data=italy, fill=NA)+
  ggtitle("Shapefile")+
  coord_sf(xlim = c(min_lon_x - error, max_lon_x + error), 
           ylim = c(min_lat_y - error, max_lat_y + error))

fig2<-ggplot()+
  geom_stars(data=dem_italy)+
  geom_sf(data=italy, fill=NA)+
  scale_fill_continuous(na.value=NA)+
  ggtitle("Raster")+
    coord_sf(xlim = c(min_lon_x - error, max_lon_x + error), 
           ylim = c(min_lat_y - error, max_lat_y + error))

Turning a Raster into a Vector

This fundamentally creates a grid

library(ggpubr)
ggarrange(fig1, fig2)

Turning a Raster into a Vector

The function st_rasterize is what allows us to convert a vector layer to a raster

This is for example the types of vectors that we can produce from rasters

Extracting Raster Values

It may be necessary to “extract” raster values according to overlapping vector features, such as points, lines or polygons.

Extracting Raster Values

It may be necessary to “extract” raster values according to overlapping vector features, such as points, lines or polygons.

Extracting Raster Values

It may be necessary to “extract” raster values according to overlapping vector features, such as points, lines or polygons.

Extracting Raster Values

It may be necessary to “extract” raster values according to overlapping vector features, such as points, lines or polygons.

# Extracting Raster Values

When extracting values to polygons, it is common to summarize the values that the geometry intersects with, using mean

Extracting to points is simpler, because each geometry corresponds to one pixel, so there is nothing to summarize.

Extracting to points can be accompished with the st_extract function

Extracting Raster Values

We can for example aggregate and take the average of all rasters at the polygon level

Loading counties in Italy

library(sf)
#Step1: Reducing the resolution: Higher dx/day values mean lower resolution
elev<-st_warp(all, st_as_stars(st_bbox(all), dx = 0.2, dy = 0.2))
#Step2: Renaming the raster
names(elev)<-"elev"
#Step3:Reading the shapefile
ita2 <- st_read(dsn="./data/gadm41_ITA_shp/gadm41_ITA_2.shp", quiet = TRUE)
#Step4:Keeping only the relevant variable
ita2<-subset(ita2, select = c(NAME_2))
#Step5: Simplify polygons
ita2<-st_simplify(ita2,  dTolerance = 0.005)
#Step6: Making sure that the raster and shapefile have the same crs
st_crs(ita2)<-st_crs(elev)
#Step7: Cropping the raster to the shapfile
dem_italy<-st_crop(elev, ita2)
#Step8: Calculate the mean by polygon
elev_by_county <- raster::aggregate(dem_italy, ita2, mean, na.rm=TRUE)
#Step9: Adding the raster by polygon
poly_elev <- mutate(ita2, avg_elev = elev_by_county$elev)

Extracting Raster Values

We can for example aggregate and take the average of all rasters at the polygon level

Loading counties in Italy

fig1<-ggplot()+
  geom_sf(data=ita2, fill=NA)+
  scale_fill_continuous(na.value=NA)+
  ggtitle("Polygons")+
    coord_sf(xlim = c(min_lon_x - error, max_lon_x + error), 
           ylim = c(min_lat_y - error, max_lat_y + error))

fig2<-ggplot()+
  geom_stars(data=dem_italy)+
  geom_sf(data=ita2, fill=NA)+
  scale_fill_continuous(na.value=NA)+
  ggtitle("Raster")+
    coord_sf(xlim = c(min_lon_x - error, max_lon_x + error), 
           ylim = c(min_lat_y - error, max_lat_y + error))

fig3<-ggplot()+
  geom_sf(data = poly_elev, mapping = aes(fill = avg_elev))+
  scale_fill_continuous(na.value=NA)+
  ggtitle("Zonal Stats")+
    coord_sf(xlim = c(min_lon_x - error, max_lon_x + error), 
           ylim = c(min_lat_y - error, max_lat_y + error))

ggarrange(fig1, fig2, fig3, ncol = 2, nrow = 2)

Extracting Raster Values

Rasterizing Polygons

We might be interested in turning the shapefiles that we just create into a raster.

We can achieve that by using st_rasterize

poly_elev_raster<-st_rasterize(poly_elev)

ggplot()+
  geom_stars(data=poly_elev_raster)+
  scale_fill_continuous(na.value=NA)

Rasterizing Polygons

We might be interested in turning the shapefiles that we just create into a raster.

We can achieve that by using st_rasterize

Rasterizing point counts

We might be interested in counting points by grid.

library(sf)
library(ggplot2)
#Step1: Reading Spain borders
esp0 <- st_read(dsn="./data/gadm41_ESP_shp/gadm41_ESP_0.shp", quiet = TRUE)
#Step2: Reading wildfires
fires <- read.csv('./data/fires_2021/viirs-snpp_2021_Spain.csv')
#Step3: Turning the acq_date variable into a date variable
fires$acq_date<-as.Date(fires$acq_date)
#Step4: Make sf object from points
fires<-st_as_sf(fires, coords = c("longitude", "latitude"))
#Step5: Adding a CRS to the shapefile
fires <- st_set_crs(fires, st_crs(esp0))  # 4326 is the EPSG code for WGS84
#Step6: Extract month and year information
fires$month <- format(fires$acq_date, "%m")
fires$year <- format(fires$acq_date, "%Y")
#Step7: Simplifying our data
fires<-subset(fires, type==0)
#Step8: Reading only August
fires_aug<-subset(fires, acq_date>"2021-07-31" & acq_date<"2021-09-01")

Rasterizing point counts

We might be interested in counting points by grid.

#Step10: Creating a grid around the country shapefile
r <- terra::rast(esp0, res=(0.5 * 0.5))
v <- terra::as.polygons(r, extent=FALSE)
#Step11: Turning the spatraster into an sf object
grid_sf<- sf::st_as_sf(v)
#Step2: Simplify shapefiles
esp0b<-st_simplify(esp0, dTolerance=2000)
#Step12: Simplifying country polygon
#esp0 = st_cast(esp0,"POLYGON")
#Step13: Calculating polygon area
esp0b$area<-st_area(esp0b)
#Step14: Turning the units into square km
esp0b$area2<-units::set_units(esp0b$area, km^2)
#Step15: Only choosing the largest polygon
esp0c<-subset(esp0b, area2==max(esp0b$area2))
#Step16: Making a box object around the polygon
bbox <- st_bbox(esp0c)
#Step17: Extracting min and max latitudes and longitudes
min_lat_y <- bbox["ymin"]
max_lat_y <- bbox["ymax"]
min_lon_x <- bbox["xmin"]
max_lon_x <- bbox["xmax"]
#Step18: Interact grid with country shape files
xsubset<-st_intersection(esp0, grid_sf)
#Step19: Mapping the cropped grid
fig1<-ggplot()+
  geom_sf(data=xsubset, fill=NA)+
  geom_sf(data=esp0, fill=NA)+
  scale_fill_continuous(na.value=NA)+
  ggtitle("Polygons")+
    coord_sf(xlim = c(min_lon_x - error, max_lon_x + error), 
           ylim = c(min_lat_y - error, max_lat_y + error))

Rasterizing point counts

#Step20: Intersect points with grid
x <- st_intersects(fires_aug, xsubset, sparse = FALSE)
#Step21: Adding the points by grid
x <- colSums(x)
#Step22: Adding the sum the grid
xsubset$no_fires<-x
#Step23: Mapping the grid
ggplot() + 
  geom_sf(data=xsubset, aes(fill=no_fires)) + 
  theme_void() +
    coord_sf(xlim = c(min_lon_x - error, max_lon_x + error), 
           ylim = c(min_lat_y - error, max_lat_y + error))

Rasterizing point counts

#Step24: Mapping the grid
ggplot() + 
  geom_sf(data=xsubset, aes(fill=log(no_fires+1))) + 
  theme_void() +
    coord_sf(xlim = c(min_lon_x - error, max_lon_x + error), 
           ylim = c(min_lat_y - error, max_lat_y + error))